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ABSTRACT 

Analytic models based on spherical and ellipsoidal gravitational collapse have been 
used to derive the mass functions of dark matter halos and their progenitors (the 
conditional mass function). The ellipsoidal model generally provides a better match 
to simulation results, but there has been no simple analytic expression in this model 
for the conditional mass function that is accurate for small time steps, a limit that 
is important for generating halo merger trees and computing halo merger rates. We 
remedy the situation by deriving accurate analytic formulae for the first-crossing dis- 
tribution, the conditional mass function, and the halo merger rate in the ellipsoidal 
collapse model in the limit of small look-back times. We show that our formulae pro- 
vide a closer match to the Millennium simulation results than those in the spherical 
collapse model and the ellipsoidal model of Sheth & Tormen (2002). 
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1 INTRODUCTION 



, iPress fc Scliecht"eil (119741 ) presented an analytical expression 
• for the unconditional mass function of dark matter halos at 
redshift z, n{M,z), based on the spherical collapse model. 
. This function is closely related to the first-crossing distri- 
' bution of r andom walks with a barrier in the excursion set 
\ framework l|Bond et al.lll99ll ). In this framework, the linear 
■ over-density computed at a given point in the Lagrangian 
' space fluctuates as a Markovian process when smoothed on 
, successively smaller scales, and a dark matter halo is iden- 
tified at the point when the random walk of the linear over- 
density crosses a critical value, or a barrier B{M, z). In the 
spherical collapse model, this barrier depends only on time 
and is independent of mass: B = 5c/D{z), where Sc ~ 1.68 
and D{z) is the standard linear growth factor. 

Although analytically simple, the spherical collapse 
model has been found to over-predict the abundance of 
s mall halos and under-predict that of massive ones (e.g. 



.Lacev fc Cold 1 19941: [Celb fc Bertschingerl 1 1994 iTormenI 
Il998l : ISheth fc Tormenlll999l ). The reason is mainly that 



halo collapses are generally triaxial rather than spherical 
(e.g. ■ [Poroshk cvich 1970; Bardcon ct al. 1986). Based on 
iBond fc MversI l| 19961 ), ISheth et all (|200ll ) use the ellipsoidal 
collapse model and obtain fitting functions that provide a 
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closer match to the unconditional halo mass function in N- 
body simulations. Unlike the spherical collapse model in 
which the condition for the virialization of a dark matter 
halo is solely determined by the linear over-density on the 
scale of the halo mass, the virialization condition in the el- 
lipsoidal collapse model also depends on halo ellipticity and 
prolateness. By assuming that a dark matter halo b e comes 
virialized when its third axis collapses, ISheth et al.l (|200ll ) 
find a new criterion for the virialization of dark matter ha- 
los, which involves all three parameters. They further sim- 
plify the virialization condition by fixing the ellipticity and 
the prolateness at their most likely values for a given over- 
density, and obtain a fitting formula for the barrier B{M, z) 
that is mass- dependent, in contrast to the constant B{z) 
of the spherical collapse model. A mass- dependent barrier 
is commonly referred to as a moving barrier, and it is this 
mass-dependence that suppresses the abundance of small 
halos while increasing that of massive ones in the ellipsoidal 
collapse model. Physically, this is because a smaller halo typ- 
ically has a larger ellipticity and therefore a longer collapsing 
time. 

The relationship between the unconditional mass func- 
tion and the first crossing distribution associated with 
barrier-crossing random walks has been exte nded to ob- 
tain the conditional ma ss function of halos l|Bond et al.l 
ll99ll : lLacev fc Colel 'l993Y In this so-called extended Press- 
Schechter (EPS) model, the conditional mass function 
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dN{Mi, zi\Mo, zo)/dAli gives the average number of pro- 
genitor halos (of mass Mi at redshift zi) per unit mass 
associated with a descendant halo of mass Mo at redshift 
zo (21 > zo). Once determined, it can be used to generate 
merger trees of halos for many applications (e.g., galaxy for- 
mation, growth of the central black hole, reionization) using 
Monte Carlo simulations. 

The conditional mass function has a simple ana- 
lytic form in the con stant barrier spherical collapse model 
l|Lacev fc Cold Il993t ). For a moving barrier (such as the 
ellipsoidal collapse model), however, exact analytic forms 
have been found only for t he special case of a linear bar- 
rier (jSheth fc TormenlbOOl ST02 hereafter); while the same 
authors have proposed a Taylor-series-like approximation 
for a general moving barrier. We find that none of these 
formulae work well for zi — zq ^ 1, which is important 
for generating a ccurate merger trees in most Monte Carlo 



meth o ds (e.g. , Lacev fc Cold 19931: Kauffmann fc Whitd 



Cole et al 



I993I: ISomerviUe fc KolattI Il999l : ISheth fc LemsonI 1 19991 : 



2OOOI) ating halo merger rates to 

the co nditional mass function (§ 2 below). iFakhouri fc Mai 
(|2007l ) compare the hal o merger rates determ ined from the 
Millennium simulation (|Springel et al.|[2005l ) with the pre- 
diction of the spherical collapse EPS model, finding the 
latter to overpredict the major merger rates by up to a 
factor of ~ 2 and underpredict the minor merger rates 
by up to a factor of ~ 5. Recently, various other fitting 
forms for the conditional mass function have been proposed 
that are calibrat ed to the results from particular A'^-bod y 
simulations, e.g. 



iNeistein fc Dekel 



Cole et all (j2007l ): IParkinson et al.1 (|2007h : 



(|2007l l. An alternative way of deriving 
the conditional mass function that does not rely explicitly 
on fitting to si mulations is to so lve the integral equation 
proposed by Zh ang fc Huil (|2006[ l (ZH06 hereafter). This 
method, which is based on the conservation of probability in 
the excursion set formalism, is accurate but computationally 
expensive. 

In §2, we derive almost exact analytic forms for the first 
crossing distribution, the conditional mass function, and the 
halo merger rate in the ellipsoidal collapse model in the 
limit of small look-back times; a limit where earlier work 
breaks down. Our method is based on ZH06, but our re- 
sults are expressed in simple analytic forms. We compare the 
predictions of this improved ellipsoidal collapse model with 
those from the spherical collapse model and the Millennium 
simulation in 331 We assume the cosmological parameters 
used in the Millennium simulation: flm = 0.25, ~ 0.045, 
h = 0.73, = 0.75, n = 1, erg = 0.9. 



2 IMPROVED ELLIPSOIDAL COLLAPSE 
MODEL FOR SMALL TIME STEPS 

We use dN{M\, zi\Mo, Zo) /dMi to denote the conditional 
mass function of dark matter halos, defined in §1. We use 
Bm{Mi, M2, z)dM-idM2 to denote the merger rate of ha- 
los, which is defined to be the number of mergers between 
halos of mass (Mi, Mi -I- dMi) and (M2,M2 -I- dM2) per 
unit volume and unit redshift at redshift z. If halo mergers 
are assumed to be binary and mass conserving such that 
the consequence of each merger is the formation of a de- 
scendant halo of mass Mo = Mi -|- M2, then we can relate 



the merger rate to the (number-weigh ted) conditional mass 
function through a simple relation (see 'Fakhouri fc Mall2007l 
for a detailed discussion, or lShethfc Pitman., 199?! for an ear- 
lier discussion): 



BMiMi, 
1 d 



M2, z)dMidM2 = n(Mi + M2, z)d{Mi + M2) 
" Az dM^ ^(^'^i' ^ + ^-^l^'^i + M2,z)dMi (1) 

where n(M, z) is the (unconditional) halo mass function at 
redshift z, and Az is assumed to be small. 

To study how many progenitors at zi are associated 
with a descendant halo of mass Mo at zo , we let the random 
walk of the linear over-density start from the scale of the de- 
scendant halo S'(Mo) with an over-density of B[S{Mo), zo], 
where S{M) — a^{M) is the variance of the linear density 
field smoothed with a window function containing mass M. 
A progenitor of mass Mi is then identified once the ran- 
dom walk crosses ^[^(Mi),^!] on the scale of S'(Mi). The 
conditional mass function can then be written as 

Ml-4^N{Ml,z^\Mo,zo)dM^ = Mof{AS)dAS (2) 
dMi 

where AS = S'(Mi) -S'(Mo), and f{AS) is the first-crossing 
distribution of random walks with a barrier of the form 
b{AS) = B[S{Mi),zi] - B[S{Mo), zo]. 

As mentioned in § 1, analytic solutions for /(AS) have 
only been found when the barrier is a constant or a linear 
function of AS. Ifere we propose an almost exact analytic 
solution for the first-crossing distribution in the limit of zi — 
Zo <^ 1. We use the fact that the barrier b{AS) is a weakly 
nonlinear function of AS to approximate it by 

b{AS) = bo + biAS + b2{ASf , (3) 

where b2{AS)'^ is assumed to be subdominant in compari- 
son with the other two term^. The first-crossing distribu- 
tion can then be written as /(AS, bo, 61, 62) • When 62 = 0, 
f{AS,bo,bi,0) has the analytic form (see ST02): 



/(AS,6o,&i,0) 



bo 



: exp 



(bo + biASf 
2AS 



(4) 



ASV2^^ 

For nonzero b2, we approximate f{AS,bo,bi,b2) as 

/(AS,bo,bi,b2) ~ /(A5,bo,bi,0) (5) 
+ b2 X dbJ{AS,bo,bi,b2)\^^^g 

and derive the second term on the right side of eq. ^ next. 

For notational simplicity, we denote AS as s and 
/(s, bo, bi, 62) as /(s). The first-crossing /(s) in general sat- 
isfies the integral equation (see ZH06 for the derivation): 



f{s)^gi{s)+ ds'f{s')g2{s,s') 
Jo 

where 



(6) 



5i(s) = 
92(3, s') = 



b_{s)_^dh 
s ds 



^db b{s)-bis') 
ds s — s' 



Polb{s),s] (7) 
Po [b{s)-b{s'),.s^s'] 



^ This is obviously not true when AS is very large. But we find 
that in practice, the relevant range of AS is rarely large enough 
to invalidate the assumption. 
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PoiS,s) 



: exp 



"2s 



From eq. (|6} and after some algebra, we find 

ds'db2f{s')g2is,s')\t2^o , 



(9) 

s')Po [bi{s — s'),s — s'] 



(10) 



5f>2/(s)|i>2=o = + / ds'db2f{s')g2is,s')\b2^o, (8) 
Jo 

where 

g(s) = -bl X Po(bo + bis,s) 

+ [ ds'f{s',bo,bi,0){bls'^ - 
Jo 

g2{s, s')j62=o = 6i X Po [bi{s - s'), i 

The two terms in g{s) come fi'om (Jbjgri and c^tj^a respec- 
tively, and a number of terms have been cancelled out us- 
ing the relation J^" ds' f{s' ,bo,bi,0)Po[bi{s — s'),s — s'] — 
Po{{bo + bis, s) that follows from eq. (|6]) for 62 = 0. 

The complicated form of g(s) in eq. ((9)1 makes it diffi- 
cult to solve eq. However, we are interested in the limit 
of small look-back times for the conditional mass function, 
which corresponds to a small barrier difference b{s). We can 
thus neglect the terms of O(feo) a-nd simplify g ao 



9is) 



^0 /,2 o^ 

— (OiS — 2) exp 



bjs 
~2~ 



We then solve for db2.f{s) by combining eqs. 
and using the Laplace transform: 



db2f[s)\b2=o = -^exp I — 



1 + 



bi^/s 

r(3/2) 



(11) 

(fTU)) . and 

(12) 



Substituting this expression back into eq. ^ , we obtain the 
first-crossing distribution for a barrier of the form b{AS) = 
bo+biAS + biiASf: 



/(A5) = 



60 



ASV2tvAS 
bob2 

— A — ^^P 



exp 



(bo + biASf 
2AS 



(13) 









r(3/2) 



+ 0{bl) 



This is our main result. We note that this equation reduces 
to the analytic expressions of eq. in the spherical collapse 
model (i.e. a constant barrier with bi — b2 — 0) and the 
ellipsoidal model with a linear barrier (i.e. 62 = 0). The 
second term on the right hand side of eq. p3p is a new 
term arising from the quadratic part of b{AS). This term is 
absent in the Taylor series approximation proposed in eq. (7) 
of ST02. The latter is obtained by replacing the barrier B 
and the variance S in their unconditional mass function in 
the ellipsoidal model with B{S{Mi), zi) - B{S{Mo), zo) and 
S{Mi) — S{Mo), respectively; that is, they assume that the 
unconditional and conditional mass functions have the same 
form. This assumption holds exactly in the spherical case 
but is not so for ellipsoidal collapse. As we illustrate in Fig. 1 

^ Note that eg. lllll l is derived by neglecting the first term on the 
right side of eq. (|9]l and approximating f{s',bo,bi,0) (eq.l|4]l) as 
o"^ X exp(— feJs/2). These approximations introduce errors 



bo/V2^ 

of order bg to g{s ^ b^) and of order 60 to g(s 
the error on db2f(s) is of order bg when s 2> 
60 when s 
smaU. 



feg). Similarly, 



bg, and of order 

6g. We find that this error is negligible when bo is 



below, their expression for the conditional mass function 
becomes inaccurate for small look-back times. 

We can now use eq. ([2)| to convert our first-crossing 
distribution in eq. p3|l into the conditional mass function. 
We first need to specify the coefficients bo,bi, and 62 for the 
barrier b{AS) = B\S(Mi) zi] - B\S(Mo)[zo ]. To do so, we 
use the barrier shape from lSheth et al] (|200ll ) and ST02 that 
has been shown to provide close fits for the unconditional 
mass function: 



B[SiM),z] = ^u;(z) [l + f3{^u)- 



(14) 



where a = 0.615, P = 0.485, 7 = 0.75, ly = io^{z)/S{M), 
uj{z) = 5c/D{z), Sc = 1.68, and D{z) is the linear 
growth factor. When zi — zo <^ 1, it is straightfor- 
ward to show that bo — AuAo, 

and 62 = --4A2/[27^5'^(A^a)]^''^ where Ao = 0.866(1 - 
0.133//o~'' "^), Al = 0.308^'o"''■"^ A2 = 0.0373^'o"°■"^ 
uo — Lu'^ (zo) / S (Mo) , and Alo — 0^(21) —uj{zo). The resulting 
conditional mass function in our ellipsoidal collapse model 
for zi — zq <^ 1 is: 



^iV(Mi,2i|Mo,2o) = -^N{Mi,zi\Mo,zo) 



(15) 

sph 



Ao exp — 



AlS 



1 + ^25" 



3/2 



1 + 



r(3/2) 



where dN/dMi^pi, = (Mo /Mi ) {dSi/dMi ) {Alo /ASV2-kAS) 
is the standard spherical model result, Ao,Ai, and A2 are 
related to vo defined above, and S — AS/ S{Mo). We find 
that neglecting the A2 term in eq. p5|l (or 62 in eq. (|13|) ) 
leads to a systematic error of > 20% for Mi ~ O.OlMo. 
It is also worth noting that since the barriers at different 
redshifts can intersect in the moving barrier model (see 
appendix A of ST02 for more details), Ao in eq. p5p can 
be negative. This occurs very rarely, however, since Ao < 
only when S{A'Io) ^ 30uj^{zo), i.e. when the descendant 
mass is much smaller than the typical halo mass at zo- 

We can now use eqs. ([T} and p5p to write down an ana- 
lytic e xpression for the halo m erger rate. Following the nota- 
tion of lFakhouri fc Mai l|2007l 'l. we express the merger rate in 
terms of the total descendant mass Mo = Mi -I- M2 and the 
mass ratio of the two progenitors ^ = M2/M1 (assuming ^ ^ 
1), and use B{Mo,^, z)dMod^ = Bm{Mi, M2, z)dMidM2 to 
relate the two rates. By treating the average merger rate per 
descendant halo, B{Mo,^, z)/n{Mo, z), as a single physical 
quantity, we find 



g(Mo,C,2) ^ B{Mo,^,z) 
n{Mo,z) 



n(Mo, 



(16) 



Ao exp 



A2S: 



sph 
q3/2 



AiS, 



1/2 



= AS,/S{Mo), = 5'(M0-S'(Mo) 
from the spherical collapse model is 



B{Mo,i,z) 



n{Mo, 



duj 



Mn" 



dS{M,] 



r(3/2) 

and the prediction 

:(17) 



sph 



dz{l + O^PU dM, AS^V2^vAS, 



where Mi can be either of the progenitors Mi or M2. We 
recall that the conditional mass function and the merger 
rate in the EPS model is not symmetric with respect to the 
two progenitor masses. This remains an unsolved problem. 
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0.01 0.1 1 10-3 0.01 0.1 

Figure 1. The conditional mass functions for the progenitor ha- 
los of a descendant halo of mass Mo = 10^^ Mq at z = 0. The 
four panels are for four look-back times: Az = 0.1,0.03,0.01, 
and 0.003. Eq. I I15I I of this paper (red solid) agrees closely with 
the exact solution from the method of ZH06 (black dotted), 
while the approximation based on eq. (7) of ST02 (blue dashed) 
overpredicts the number of progenitors for small look-back time 
{Az < 0.03). 



0.083, 02 = -2.01, as = 0.409, and at = 0.371. Their 
Fig. 15 illustrates the large discrepancy between eq. ([18} 
and the prediction of the standard spherical EPS model. 

Figs. [2] (for option A) and [3] (option B) show the ratio 
between the halo merger rate B/n of our ellipsoidal col- 
lapse model (eq. (|16p ) and that of the Millennium simu- 
lation (ea. (|18p ) for three descendant halo masses and four 
redshifts. The spherical collapse model (eq. (117^ 1 is also 
shown for comparison. The minimum halo mass is chosen 
to be 2 X 10^° as set by the halo mass resolution in 
the Millennium simulation. Comparison of the two figures 
shows that the two choices of Adi give similar results for ma- 
jor mergers but yield very different predictions for ^ ^ 1, 
where option A (Fig. [2]) agrees better with the Millennium 
than option B (Fig. We note that the two options pre- 
dict different power- law dependencies on ^ at ^ <C 1: for 
S{M) oc M~^, option A gives B/n oc ^^/^-^^ but option 
B gives B/n oc which is independent of the density 

variance on the scale of the smaller progenitor mass. It is 
also interesting to note that option A is implicitly used in 
some Monte Carlo methods; for example, ,Cole e t al. (200^ 
select the mass of the first progenitor from the lower half of 
the conditional mass function (i.e. Mi < Mo/2). 

According to Fig. [2l the discrepancy between our ellip- 
soidal collapse model and the Millennium simulation is typi- 
cally 20 ~ 30%, but can reach up to about 80% when ^ ~ 1. 
On the other hand, the relative difference between the spher- 
ical collapse model and Millennium is typically 40% ~ 60%, 
and reaches up to 120% for the major mergers. Therefore 
in almost every case, the new ellipsoidal collapse model im- 
proves the agreement with the Millennium simulation. 



Below we simply show the results for b o th cho ices. Accord- 
ing to the notation of iFakhouri fc Mai (|2007l ). we call the 
merger rate "option A" when Mi is assigned to the smaller 
progenitor M2, and "option B" when Mi — Mi. 



3 NUMERICAL RESULTS AND COMPARISON 
WITH THE MILLENNIUM SIMULATION 

In Fig. [1] we illustrate the accuracy of eq. (|15|l by comparing 
it with the exact solution from ZH06 and the analytic ap- 
proximation based on eq. (7) of ST02. The figure shows the 
conditional mass functions at four different look-back times 
{Az = 0.1,0.03,0.01,0.003) for a descendant halo of mass 
10" M0 at redshift zero. The approximation of ST02 is seen 
to overpredict the number of lower mass progenitors by up 
to a factor of 2 to 10 for Az < 0.03, whereas eq. (|15|) of this 
paper is accurate when Az is as small as 0.003. 

For the halo merger rate, we compare the predictions of 
our ellipsoidal mod el with the Millennium simulation merger 
rate determined bv IFakhouri fc Mai ()2007[ ) using the "stitch- 
ing" method. They find that the halo merger rate in the 
simulation converges well when the look-back time Az ap- 
proaches zero and can be described by a simple universal 
fitting formula: 



B{Mo,^,z) 



\M J 



n{Mo,z 
where M = 1.2 x lO^^M©, | 



Car> 
exp 



0.098, A 



duj{z) 



dz 

0.0289, ai 



(18) 



4 SUMMARY 

We have derived new analytic formulae for the first- 
crossing distribution (eq.((T3}), the conditional mass function 
(ea. (|15|l ). and the merger rate of dark matter halos (eq. (|16p ) 
in the ellipsoidal collapse model in the limit of small look- 
back times. Our method is based on solving the first-crossing 
distribution of random walks with a moving barrier using 
the exact integral equation of ZH06. This method results in 
extra terms in eqs. (|13|l . (|15p . and (|16p that are absent in 
the spherical collapse model and different from those in the 
ellipsoidal model of ST02. Fig. 1 illustrates how these terms 
correct the discrepancies of ST02 in the conditional mass 
function for small look-back times Az < 0.03. 

Eq. ((T} relates the conditional mass function at small 
Az to the halo merger rate. The halo merger rate of our 
ellipsoidal collapse model gen erally agrees bet t er wit h the 
Millennium result reported in IFakhouri fc Mai l|2007l ) than 
that of the spherical collapse model (red vs blue curves 
in Figs. [2] and [3]). The discrepancy between our ellipsoidal 
collapse model and the Millennium simulation is typically 
20% ~ 30%, which is about a factor of two smaller than 
that of the spherical collapse model. A comparison between 
Figs. (21 and 13 shows a better agreement between model and 
simulation when the progenitor mass Mi in the analytic for- 
mulae is assigned to be the smaller progenitor (option A). 

A number of factors not considered in this paper nor 
in earlier EPS work can contribute to the remaining 20-30% 
discrepancy between the model and simulation. These in- 
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Figure 2. Comparison of the halo merger rates from our eUip- 
soidal model predictions (ea. lll6ll7ll ) vs. the Millennium simula- 
tion f eg. Hist ). The four panels show the merger rates as a function 
of the progenitor mass ratio ^ = M2/M1 at 2 = 0, 1,2,3. Within 
each panel, the red solid curves show the ratio of our ellipsoidal 
collapse model prediction (eq. JTHJ) to the Millennium result; the 
blue dashed curves show the ratio of the standard spherical col- 
lapse model (eg.l llTIl ) to the Millennium result. For each colour, 
the set of three curves show three descendant halo masses (from 
bottom up): lO^^, lO^^, and IO^^Mq. The progenitor mass Mi 
in eqs. I|17|l and ||16|I is chosen to be the less massive M2 (option 
A), which we find to match the simulation better than option B. 



elude the statistical importance of non-binary mergers, dif- 
fuse accretion and tidal stripping of halo mass (such that 
Afo 7^ Ml + M2), and the impact of a halo's environment on 
merger rates that can lead to non -Markovian processes i n 
the excursion set model {e.g. ST02. iNeistein fc Dekeill2007l ). 

Our conditional mass function in ea. (|15|l can be eas- 
ily incorporated into Monte Carlo simulations to study 
halo merger histories over a large look-back time. This is 
done in a companion paper, in which we compare sev- 
eral existing Monte Carlo algor i thms (e.g. . iLacev fc Colel 
I993I: iKauffmann fc Wh"it3l 19931 : ISomerville fc Kolattll 19991 : 



Cole::kalJbQd) that ~~b ;sed on the spheric; n^: 



lapse model, and propose a more accurate method using 
ea. (|15p . Several groups have recently proposed accurately 
parameterized forms of the conditional mass function based 
on fits to the Millennium results and i ncorporated them 
into different Monte Carlo simulati ons fcol e et al.l l2007l : 
[Parkinson et al.l l2007l : INeistein fc De kel 2007) . As we will 
show in the companion paper (also see lMoreno fc Shethl2007l 
for a different method), a similar level of accuracy can be 
achieved using eq. (|15p of this paper and our Monte Carlo 
method without a priori knowledge of N-body simulations. 

We thank James Bullock, Joanne Cohn, Lam Hui, Ravi 
Sheth, Martin White, and Simon White for useful discus- 
sions. This work is supported in part by NSF grant AST 



Figure 3. Same as Fig. [21 except that the progenitor mass Mi in 
both cq. JTtJ and I II6I I is chosen to be the more massive progenitor 
Ml (option B), which is not our preferred option. 



0407351. The Millennium Simulation databases used in this 
paper and the web application providing online access to 
them were constructed as part of the activities of the Ger- 
man Astrophysical Virtual Observatory. 
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